How the storage effect and the number of temporal niches affect biodiversity in stochastic and seasonal environments

The main metric for the persistence of a two-species community is T, the mean time until one of the species is lost. We will refer to T as the time to extinction (of either species).

In our two-species community the state of the system at a given time is fully characterized by the frequency of the focal species x and the state of the environment. Since we consider only the annealed dynamics where x variations are much slower than environmental variations, we characterize the initial state of the community by x alone, so T(x) is the mean time until extinction, averaged over all initial states of the environment and over all histories [38].

Implementing the diffusion approximation one finds that T satisfies [38–40], (8) where μ(x) is the mean velocity, σ2(x) is the associated variance (see Methods, Section V.A) and prime represent a derivative with respect to x. Note that in Eq (8) time is measured in elementary steps, to translate that into generations T must be divided by N.

In the Methods section we calculate the values of μ(x) and σ2(x) for four combinations of global (G) and local (L) competition with stochastic (S) and periodic (P) environmental variations.

Local-periodic: As shown in the previous section, when the competition is local dx/dt is an odd function of Δs. As explained in the Methods section, this implies that only s0 contributes to the velocity. Moreover, if the dynamic is periodic γ and δ have no effect on σ2(x) and μ(x), which are simply those obtained in the standard case of a fixed environment, (9)

Global-periodic: when the competition is global Eq (7) suggests a bias towards the coexistence point. This bias is independent of the sign of Δs and its strength is proportional to (Δs)2, which, when the diffusion approximation holds, may be approximated by γ2. In the Methods section we found an additional term in the expression for μ(x). This new term is proportional to γ2 and represents a bias towards x = 1/2. (10)

Since the variations are periodic, they do not contribute to the diffusion term σ2(x).

Local-stochastic: stochastic abundance variability increases due to the chance of the focal species, say, to pick a sequence of good or bad years. In the local stochastic case the resulting expressions (Methods, V.A) are (g is defined in Table 1), (11)

As before, the μ(x) has a term that appears to provide a bias towards x = 1/2. However, the diffusivity σ2(x) is maximal at x = 1/2 and vanishes close to the extinction point x = 0 and x = 1. The system tends to stick to the regions where its diffusion constant is small [6], and this “diffusive trapping” acts against the stabilizing effect of the μ term. These two contradicting tendencies are known to cancel each other exactly [6, 41]. As a result, the only net effect of stochastic environmental variations is an increase in the amplitude of abundance variations, which decreases extinction times [32, 38].

Global-stochastic: finally, when competition is global the (Δs)2 term facilitates coexistence. As shown in Methods section, the relevant terms are, (12)

While the diffusion term is the same as in Eq (11), the velocity has an extra piece that promotes stability, γ2x(1 − x)(1 − 2x)/2. This breaks the tie in favor of stabilization and the storage effect manifests itself [33, 38].

Eq (8), with the relevant μ(x) and σ2(x), is an inhomogeneous and linear first order differential equation for T′, so one may solve it via an integration factor and then find T by an additional integration. When this procedure is implemented numerically, as detailed in the Methods section (V.B), the results fit perfectly the outcomes of our Monte-Carlo simulations, as demonstrated in Fig 1.

Black circles are the outcomes of a Monte-Carlo simulation, and the colored lines are the theoretical predictions for the corresponding case (see legend) as obtained from numerical solutions of Eq (8) with the μ(x) and σ2(x) from Eqs (9)–(12). Details of the solution are presented in the Methods section. The hierarchy TGP > TGS > TLP > TLS is evident. In panel (a), both species have the same mean fitness (s0 = 0) and therefore all the lines are symmetric around x = 1/2. Other parameters in panel (a) are N = 600, γ = 0.25 and δ = 0.55. In panel (b) the focal species is slightly advantageous, with s0 = 0.01 (all other parameters are the same), so the maximum time to extinction appears at x < 1/2.

https://doi.org/10.1371/journal.pcbi.1009971.g001

Fig 1 reveals the hierarchy of stability properties. The local-stochastic case has the shortest time to extinction, as the role of environmental variations is purely destabilizing. In the local-periodic case this destabilizing effect of environmental variations averages out over each cycle, so its T is larger. Under global competition temporal variations facilitate coexistence, so their T is larger than in the local cases. Finally, time to extinction in the global-stochastic scenario is shorter than in its global-periodic counterpart: a sequence of bad years may cause extinction in the former, but not in the latter case.

More quantitatively, the mean time to extinction in the local-stochastic and in the global-stochastic cases, and in addition T in fixed environment (which, as explained, is equivalent to the local-periodic case), were calculated by [38]. Let us quote some of their result, and contrast them with the new results obtained in the global-periodic case.

When s0 = 0, or otherwise in the weak selection regime where the effect of s0 is negligible, the maximum value of T (max over all values of initial state x) is,

To understand these different dependencies on N, let us neglect the effect of demographic stochasticity and replace it with an absorbing threshold at z values that correspond to a single individual, z ≈ ±ln N. In the local-stochastic case with s0 = 0 the abundance performs an unbiased random walk along the logit (z) axis, therefore exit times scale with ln2 N [40]. When competition is global and the environment is stochastic the μ(x) term supports an attractive fixed point in x = 1/2, and extinction occurs through improbable sequences of bad years. During such a sequence abundance decreases exponentially, so the number of bad years required to cross the one individual threshold is proportional to ln N. Since the chance of such a sequence decreases exponentially with its length, the mean time to extinction scales like a power-law in N.

Finally, in the global-periodic case (that was not discussed by [38]) there are no such rare sequences of bad years. Therefore, extinction may take place only due to demographic stochasticity (note that the σ2(x) term in Eq (10) is the same as in Eq (9), reflecting only demographic variations). This requires a highly improbable sequence of death events that may take a population with order N individuals to extinction, and such a sequence is exponentially rare in N. Therefore, one expects T in the global-periodic scenario to grow exponentially with N. These four behaviors: ln2 N, linear, power-law, and exponential, are shown in Fig 2.

Filled red circles were obtained from Monte-Carlo simulations, Blue circles are theoretical predictions based on numerical integration of Eq (8) with the μ(x) and σ2(x) from Eqs (9)–(12), and the dashed black line is a linear fit to these blue circles, presented to guide the eye. As expected, in the global-periodic case the mean time to extinction grows exponentially with N, in the global-periodic the grows satisfies a power-law, the local-periodic case behaves like the neutral model (T is linear in N), whereas the local-stochastic dynamics yields log2 N growth. In all cases s0 = 0 (so the T shown here is the mean time to extinction starting from x = 1/2), δ = 0.2 and γ = 0.4.

https://doi.org/10.1371/journal.pcbi.1009971.g002

When one cannot neglect s0 the situation is quantitatively different, but the qualitative hierarchy is preserved [38]). In that case under local competition T ∼ ln N (in both periodic and stochastic cases). When competition is global, as long as |s0| is not too large and μ(x) still posses an attractive fixed point at some 0 < x < 1, extinction takes place via accumulation of bad years (in the stochastic case) or random death events (in the periodic case). As a result T has the same scaling with N, although the actual time to extinction for a given set of N, γ and δ become shorter as |s0| grows.